#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <getopt.h>
#include <sys/time.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>

#include "mc-option.h"

int main(int argc, char *argv[]) {
    mcPars pars;
    mcResults results;
    struct timeval start, end;
    double runtime;
    
    getPars(argc, argv, &pars);
    
    if(1337) {
        double d1, d2, callValue, putValue, delta;

        d1 = (log(pars.S0 / pars.strikePrice) + (pars.r + pars.sigma*pars.sigma * 0.5) * pars.endTime)
                / (pars.sigma * sqrt(pars.endTime));
        d2 = d1 - pars.sigma * sqrt(pars.endTime);
        callValue = pars.S0 * gsl_cdf_ugaussian_P(d1) - pars.strikePrice * exp(-pars.r * pars.endTime) * gsl_cdf_ugaussian_P(d2);
        putValue = pars.strikePrice * exp(-pars.r * pars.endTime) * gsl_cdf_ugaussian_P(-d2) - pars.S0 * gsl_cdf_ugaussian_P(-d1);
        delta = gsl_cdf_ugaussian_P(d1);
        if(pars.isPut) {
            printf("Black-scholes value for european put option: %f, delta %f\n", putValue, delta);
        } else {
            printf("Black-scholes value for european call option: %f, delta %f\n", callValue, delta);
        }
        
        double N = pars.endTime * pars.checksPerYear;
        d1 = sqrt((6.*N*N) / ((N+1.) * (2.*N+1))) * (log(pars.S0 / pars.strikePrice) 
            + pars.r * (N+1.)/(2.*N) * pars.endTime + (N+1.)*(N+2.)/(12.*N*N) * pars.sigma * pars.sigma * pars.endTime)
                / (pars.sigma * sqrt(pars.endTime));
        d2 = d1 - sqrt((N+1.) * (2.*N+1)/(6.*N*N))*pars.sigma * sqrt(pars.endTime);
        callValue = pars.S0 * exp((1.-N)/(2.*N) * pars.r * pars.endTime + (1.-N*N)/(12.*N*N)*pars.sigma * pars.sigma *pars.endTime)
        * gsl_cdf_ugaussian_P(d1) - pars.strikePrice * exp(-pars.r * pars.endTime) * gsl_cdf_ugaussian_P(d2);
        pars.asianValue = callValue * exp(pars.r * pars.endTime);
        printf("Black-scholes value for geom. asian option: %f\n", callValue);
    }   
    
    if(pars.bumpEps > 0.0) {
		int k;
		double S0 = pars.S0;
		double x;
		double mean=0.0, M2=0.0, variance, dVar, stdDev;
		for(k=1; k<=200; k++) {
			mcOption(&pars, &results);
			x = results.optionValue;
			pars.S0 = S0 + pars.bumpEps;
			if(pars.bumpRenewSeed)
				pars.seed++;
			mcOption(&pars, &results);
			x = (results.optionValue - x) / pars.bumpEps;
			dVar = x-mean;
			mean += dVar/k;
			M2 += dVar * (x-mean);
			
			pars.S0 = S0;
			pars.seed++;
		}
		variance = M2/(k-1);
		stdDev = sqrt(variance);
		printf("Hedge delta\tStd dev  \tin percent  \tStd error\n");
		printf("%lf\t%e\t%lf%%\t%e\n", mean, stdDev, stdDev/mean * 100.0, sqrt(variance/k));
	} else {
		gettimeofday(&start, NULL);
		mcOption(&pars, &results);
		gettimeofday(&end, NULL);
		runtime = end.tv_sec - start.tv_sec +  (end.tv_usec - start.tv_usec) * 1.e-6;
		printf("---\nOption value\tStandard error\tHedge delta\tTime\n");
		printf("%lf\t%e\t%lf\t%lf\n", results.optionValue, results.stdErr, results.hedge, runtime);
	}
    
    
    return 0;
}

void mcOption(mcPars *pars, mcResults *r) {
    int k;
    double M2 = 0.0, mean = 0.0, dVar, variance, x;
    double confint;
    int seed = pars->seed;
    
    gsl_rng *randomNumGen = gsl_rng_alloc(gsl_rng_mt19937);
    //gsl_rng *randomNumGen = gsl_rng_alloc(gsl_rng_ranlxs2);
    gsl_rng_set(randomNumGen, seed);
    
    for(k=1; k <= 0.5 * pars->nTrials; k++) {
        if(pars->controlVariates == 0) {
            x = stockSim(pars, randomNumGen);
        } else {
            x = stockSimCtrl(pars, randomNumGen);
        }
        dVar = x-mean;
        mean += dVar/k;
        M2 += dVar * (x-mean);
        
        if(pars->doPrint && k % pars->stepsPerOutput == 0) {
            variance = M2/(k-1);
            confint = sqrt(variance / (double)k);
            printf("k: %8i\tVariance: %e\tMC stdDdev: %e\n", k, variance, confint);
        }
    }
    variance = M2/(k-1);
    r->hedge = 0.0;
    r->stdErr = sqrt(variance / (double)k);
    r->optionValue = exp(-pars->r * pars->endTime)*mean;
}


// Simulate one option up to maturity and return the value
double stockSim(mcPars *pars, gsl_rng *randomNumGen) {
    int ctr;
    double S1, payoff1, totalPayoff, rand, rand2;
    totalPayoff = 0.;

    if(pars->isAmerican == 0 && pars->isAsianA == 0 && pars->isAsianG == 0) { // We can instantly go to maturity
        rand = gsl_ran_ugaussian(randomNumGen);
        if(pars->antithetic == 1) {
            rand2 = -rand;
        } else {
            rand2 = gsl_ran_ugaussian(randomNumGen);
        }
        
        for(ctr = 0; ctr < 2; ctr++) {  //Do two paths
            S1 = pars->S0 * exp( (pars->r - 0.5*pars->sigma*pars->sigma) * pars->endTime
                            + pars->sigma * rand * sqrt(pars->endTime));
            if(pars->isDigital) {
                if(pars->isPut) {
                    payoff1 = (S1 < pars->strikePrice) ? 1.0 : 0.0;
                } else {
                    payoff1 = (S1 > pars->strikePrice) ? 1.0 : 0.0;
                }
            } else {
                if(pars->isPut) {
                    payoff1 = fmax(pars->strikePrice - S1, 0.0);
                } else {
                    payoff1 = fmax(S1 - pars->strikePrice, 0.0);
                }
            }
            if(pars->likelihoodRatio) {
				payoff1 = payoff1 * rand / (pars->S0 * pars->sigma * sqrt(pars->endTime));
			}
            rand = rand2;
            totalPayoff += payoff1;
        }
    } else { // We need to calculate the path
        int i, nsteps;
        double dt, sum1, sum2, geomAvg1, geomAvg2;
        double payoff2, S2;
        nsteps = round(pars->endTime * pars->checksPerYear);
        dt = pars->endTime / nsteps;

        S1 = pars->S0;
        S2 = pars->S0;
        sum1 = 0.;
        sum2 = 0.;
        geomAvg1 = 1.;
        geomAvg2 = 1.;
        
        for(i=0; i<nsteps; i++) {
            rand = gsl_ran_ugaussian(randomNumGen);
            if(pars->antithetic == 1) {
                rand2 = -rand;
            } else {
                rand2 = gsl_ran_ugaussian(randomNumGen);
            }
            
            S1 *= exp( (pars->r - 0.5*pars->sigma*pars->sigma) * dt
                            + pars->sigma * rand * sqrt(dt));
            S2 *= exp( (pars->r - 0.5*pars->sigma*pars->sigma) * dt
                + pars->sigma * rand2 * sqrt(dt));   
            
            if(pars->isAsianA) {
                sum1 += S1;
                sum2 += S2;
            } else if(pars->isAsianG) {
                geomAvg1 *= pow(S1, 1./nsteps);
                geomAvg2 *= pow(S2, 1./nsteps);
            }
        }

        if(pars->isAsianA) {
            payoff1 = sum1 / nsteps - pars->strikePrice;
            payoff2 = sum2 / nsteps - pars->strikePrice;
        } else if(pars->isAsianG) {
            payoff1 = geomAvg1 - pars->strikePrice;
            payoff2 = geomAvg2 - pars->strikePrice;
        }

        if(pars->isPut) {
            payoff1 *= -1.;
            payoff2 *= -1.;
        }
        payoff1 = fmax(payoff1, 0);
        payoff2 = fmax(payoff2, 0);
        totalPayoff = payoff1 + payoff2;
    }
    return 0.5 * totalPayoff;
}

// Simulate one option up to maturity and return the value
double stockSimCtrl(mcPars *pars, gsl_rng *randomNumGen) {
    int i, n;
    double S1, payoff1, totalPayoff, rand, rand2;
    static int ctr = 0;
    static double controlBeta;
    static double *x, *y;
    
    if(ctr == 0) {
        x = malloc(pars->nTrials * sizeof(double));
        y = malloc(pars->nTrials * sizeof(double));
        int i, nsteps;
        double dt, sum1, sum2, geomAvg1, geomAvg2;
        double payoff2, S2;
        nsteps = round(pars->endTime * pars->checksPerYear);
        dt = pars->endTime / nsteps;
        for(n = 0; n <= pars->nTrials * 0.5; n++) {
            // We need to calculate the path
            S1 = pars->S0;
            S2 = pars->S0;
            sum1 = 0.;
            geomAvg2 = 1.;
            
            for(i=0; i<nsteps; i++) {
                rand = gsl_ran_ugaussian(randomNumGen);
                if(pars->antithetic == 1) {
                    rand2 = -rand;
                } else if(pars->controlVariates) {
                    rand2 = rand;
                } else {
                    rand2 = gsl_ran_ugaussian(randomNumGen);
                }
                
                S1 *= exp( (pars->r - 0.5*pars->sigma*pars->sigma) * dt
                                + pars->sigma * rand * sqrt(dt));
                S2 *= exp( (pars->r - 0.5*pars->sigma*pars->sigma) * dt
                    + pars->sigma * rand2 * sqrt(dt));   
                sum1 += S1;
                geomAvg2 *= pow(S2, 1./nsteps);
            }

            payoff1 = sum1 / nsteps - pars->strikePrice;
            payoff2 = geomAvg2 - pars->strikePrice;

            if(pars->isPut) {
                payoff1 *= -1.;
                payoff2 *= -1.;
            }
            payoff1 = fmax(payoff1, 0);
            payoff2 = fmax(payoff2, 0);
            x[n] = payoff1;
            y[n] = payoff2;
        }
        //Now compute coefficient
        double vary = 0., avgy = 0., avgx = 0., covxy = 0.;
        for(i = 0; i <= pars->nTrials * 0.5; i++) {
            avgy += y[i];
            avgx += x[i];
        }
        avgy /= (double)(pars->nTrials * 0.5);
        avgx /= (double)(pars->nTrials * 0.5);
        for(i = 0; i <= pars->nTrials * 0.5; i++) {
            vary += (y[i]-avgy)*(y[i]-avgy);
        }
        for(i = 0; i <= pars->nTrials * 0.5; i++) {
            covxy += (y[i]-avgy) * (x[i]-avgx);
        }
        controlBeta = - covxy / vary;
        printf("controlbeta = %lf\n", controlBeta);        
    }
    totalPayoff = x[ctr] + controlBeta * (y[ctr] - pars->asianValue);
    ctr++;
    return totalPayoff;
}

void getPars(int argc, char *argv[], mcPars *pars) {
    pars->isPut = 0;
    pars->isAmerican = 0;
    pars->isAsianA = 0;
    pars->isAsianG = 0;
    pars->isDigital = 0;
    pars->antithetic = 0;
    pars->controlVariates = 0;
    pars->likelihoodRatio = 0;
    pars->controlBeta = 0.99;
    pars->checksPerYear = 365;
    pars->doPrint = 0;
    pars->endTime = 1.;
    pars->bumpRenewSeed = 0;
    pars->nTrials = 10000;
    pars->seed = 31337;
    pars->stepsPerOutput = 1000;
    pars->strikePrice = 99.;
    pars->S0 = 100.;
    pars->sigma = 0.2;
    pars->r = 0.06;
    pars->bumpEps = 0.00;
    
    int ch;
    while ((ch = getopt(argc, argv, "pAGDaUVvRLb:n:o:t:k:S:s:r:c:x:e:")) != -1) {
        switch(ch) {
        case 'p': pars->isPut = 1; break;
        case 'U': pars->isAmerican = 1; break;
        case 'A': pars->isAsianA = 1; break;
        case 'G': pars->isAsianG = 1; break;
        case 'D': pars->isDigital = 1; break;
        case 'L': pars->likelihoodRatio = 1; break;
        case 'a': pars->antithetic = 1; break;
        case 'V': pars->controlVariates = 1; break;
        case 'v': pars->doPrint = 1; break;
        case 'R': pars->bumpRenewSeed = 1; break;
        case 'n': pars->nTrials = (int)strtod(optarg, 0); break;
        case 'c': pars->checksPerYear = strtol(optarg, 0, 10); break;
        case 'o': pars->stepsPerOutput = strtol(optarg, 0, 10); break;
        case 'x': pars->seed = strtol(optarg, 0, 10); break;
        case 'b': pars->controlBeta = strtod(optarg, 0); break;
        case 't': pars->endTime = strtod(optarg, 0); break;
        case 'k': pars->strikePrice = strtod(optarg, 0); break;
        case 'S': pars->S0 = strtod(optarg, 0); break;
        case 's': pars->sigma = strtod(optarg, 0); break;
        case 'r': pars->r = strtod(optarg, 0); break;
        case 'e': pars->bumpEps = strtod(optarg, 0); break;
        }
    }
    
    printf("Parameters:\
            \n isPut (-p)       \t%d\
            \n isAmerican (-U)  \t%d\
            \n isAsianA (-A)    \t%d\
            \n isAsianG (-G)    \t%d\
            \n isDigital (-D)   \t%d\
            \n antithetic (-a)   \t%d\
            \n controlVariates (-V)\t%d\
            \n checksPerYear (-c)\t%d\
            \n doPrint (-v)     \t%d\
            \n nTrials (-n)     \t%d\
            \n seed (-x)        \t%d\
            \n stepsPerOutput (-o)\t%d\
            \n endTime (-t)    \t%.2e\
            \n strikePrice (-k)\t%.2e\
            \n S0 (-S)         \t%.2e\
            \n sigma (-s)      \t%.2e\
            \n r (-r)          \t%.2e\
            \n bumpEps (-e)    \t%.2e\
            \n controlBeta (-b)      \t%.2e\n\n",
            pars->isPut,
            pars->isAmerican,
            pars->isAsianA,
            pars->isAsianG,
            pars->isDigital,
            pars->antithetic,
            pars->controlVariates,
            pars->checksPerYear,
            pars->doPrint,
            pars->nTrials,
            pars->seed,
            pars->stepsPerOutput,
            pars->endTime,
            pars->strikePrice,
            pars->S0,
            pars->sigma,
            pars->r,
            pars->bumpEps,
            pars->controlBeta);
}

